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We investigate the phenomenology of freely expanding fluids, with different material properties, 
evolving through the Israel-Stewart (IS) causal viscous hydrodynamics, and compare our results 
with those obtained in the relativistic Eckart-Landau-Navier-Stokes (ELNS) acausal viscous hydro- 
dynamics. Through the analysis of scaling invariants we give a definition of thermalization time 
which can be self-consistently determined in viscous hydrodynamics. Next we construct the solu- 
tions for one-dimensional boost-invariant flows. Expansion of viscous fluids is slower than that of 
one-dimensional ideal fluids, resulting in entropy production. At late times, these flows are reason- 
^— ^ , ably well approximated by solutions obtained in ELNS hydrodynamics. Estimates of initial energy 

' densities from observed flnal values are strongly dependent on the dynamics one chooses. For the 

' same material, and the same final state, IS hydrodynamics gives the smallest initial energy den- 

sity. We also study fluctuations about these one-dimensional boost-invariant backgrounds; they are 
damped in ELNS hydrodynamics but can become sound waves in IS hydrodynamics. The difference 

Qis obvious in power spectra due to clear signals of wave-interference in IS hydrodynamics, which is 
completely absent in ELNS dynamics. 

(N 
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^ I. INTRODUCTION 

^3 . 

^ . Hydrodynamics is an effective long-distance theory of fluids close to thermal and chemical equilibrium. The hydro- 
dynamic modes are the fields which enter the theory, and can be identified from the symmetries of the microscopic 
' theory. They include the energy-momentum tensor and densities of conserved quantities such as the baryon number 
^ ■ and electric charge. The theory contains several unknown quantities, the transport coefficients and relaxation times, 
00 [ which one must obtain from microscopic computations in kinetic theory or through measurements. 

i The relativistic version of the Navier-Stokes equation was first explored by Eckart [2| and subsequently by Landau 
' 3 , who developed what we call the ELNS formalism. The relativistic equations for ideal fiuids arc widely used in the 
, contexts of heavy-ion collisions [1, 0, H and cosmology [1| . It has been suggested that the fiuid produced in heavy- 
ion collisions at RHIC is very close to ideal. Such a claim must be substantiated by a study of viscous fiuid dynamics. 
Until now there have been very few studies of non-ideal fluids in this context [T^, [ill, [Hj [Hi [13, [S [IB, [13 ■ 

ELNS theory for non-ideal fluids is known to violate causality [18,] . The problem can be traced to the linear relation 
between fluxes and thermodynamic forces which is inherent in the Chapman-Enskog method of obtaining the ELNS 
equations from kinetic theory. It was realized by Israel and Stewart [l9| that the problem with causality could be 
repaired by simply going beyond the linear relation between fluxes and forces. Their formulation of hydrodynamics 
is variously known as second-order theory, causal viscous hydrodynamics, or Israel-Stewart (IS) hydrodynamics. It 
' contains an expanded set of material parameters when compared to ELNS theory. 

There are other attempts to repair the loss of causality in Navier-Stokes theory. It was shown that the introduction 
of a phenomenologically motivated lag between the application of a thermodynamic force and the material response, 
through a memory kernel, could restore causality [201] . Such a phenomenological approach contains a smaller number 
of material parameters than IS hydrodynamics. In any case, little is known about some of the new parameters which 
appear in IS hydrodynamics. Even the relation between current correlation functions in a thermal quantum field 
theory and these quantities [2l| has not been studied comprehensively. 

In fact, once the problems of principle were resolved, further investigations of causal viscous hydrodynamics lagged 
because of an apparent paucity of applications. One expects that the main applications would be in situations where 
either the mean free path is comparable to the size of the region of interest or the relaxation time approaches the time 
scale of the phenomena of interest. It turns out that such applications are not hard to come by today. The former 
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are possibly realized in heavy- ion collisions [H, [H, [31 1 and the diffusion of neutrinos through supernovae, the latter 
in astrophysical shock waves and freezeout in relativistic reactive fluids. Interest in such systems is on the increase. 

Very little is presently known about the nature of fluid flows in IS hydrodynamics. The present paper is a s tep 
towards repairing this neglect, keeping future applications to heavy-ion physics in mind. In common with pH [l^. [isl. 

we investigate the equations for a fluid characterized completely by the energy-momentum tensor, i.e., neglecting 
net baryon number and electric charge, keeping only the shear part of the viscous stress tensor. We set up the 
equations in curvilinear coordinates appropriate to the approximate boost-invariant geometry of heavy-ion collisions, 
and reduce the tensor equations to coupled scalar equations. By examining the symmetries of the hydrodynamic 
equations we obtain laws of physical similarity. 

Since an ideal fluid has zero viscosity and zero mean free path, the ideal hydrodynamic equations are unable to 
predict their own failure. All estimates of thermalization and freezeout in ideal hydrodynamics are imposed from 
outside. Since causal viscous hydrodynamics contains the relaxation time, r^, one expects to improve upon this. We 
present a preliminary analysis. 

Next we analyze boost-invariant solutions with three sets of constitutive relations for the fluid. In each case, we 
compare the ELNS and IS descriptions of boost-invariant flow and find that the former approximates the latter at 
late times. Our results for the case of the massless Boltzmann fluid are consistent with those given in [l^, 
We also examine the propagation of fluctuations around the boost- invariant solutions. Here the differences between 
ELNS and IS descriptions are remarkable: no propagating solutions exist in ELNS dynamics, whereas IS dynamics 
gives rise to damped sound waves. 

The plan of the paper is as follows. The next section introduces the hydrodynamic equations and extracts scaling 
laws from them. This section also contains a discussion of the properties of materials that are needed in the remainder 
of the paper. Following this, we present the well-known Bjorken solution, and illustrate our method with the analysis 
of fluctuations around this solution for the ideal fluid. The next three sections examine a simple fluid, a Boltzmann 
fluid, and conformal fluids, respectively. In each case we examine the boost-invariant solution in IS hydrodynamics, 
its approach to ELNS hydrodynamics at late times, entropy production, and fluctuations around the boost-invariant 
solution. Through this analysis we build up a picture of general properties of the flow, as opposed to those which 
are speciflc to certain kinds of fluids. Wc summarize our understanding in the final section. The appendices contain 
details of the tensor analysis, the reduction of the tensor hydrodynamic equations to coupled scalar equations, and 
an analysis of transients. 



II. THE HYDRODYNAMIC EQUATIONS 



Throughout this paper we investigate the hydrodynamic equations in the limit of zero net quarks, i.e., vanishing 
baryon and charge density, since this is a good approximation to the actual situation realized in ultra-relativistic 
heavy- ion collisions. We also ignore the heat flux as in (Til . [T^ . [isl . [T3 |. References [l^, [2^ extend this to the case 
where the heat flux, thermal conductivity and baryon density are all included. It is strongly suspected that in the 
high-temperature phase of QCD, not very close to the crossover temperature Tc, the bulk viscosity is negligible. In 
view of this we investigate the equations where the bulk viscosity has been set to zero. With these simpliflcations the 
explicit form of the equations of Israel-Stewart (IS) hydrodynamics was written down in pTi. [isl [l^. 

Having set the baryon and charge densities to zero, one is left with three independent hydrodynamic variables: a 
scalar, a vector and a tensor. The scalar is the energy density, e. It is related to the pressure, p, by the equation of 
state — p — c^e, where Cs is the speed of sound. Since all three quantities in the equation of state can be written as a 
function of the temperature, T, we sometimes trade e for T. This variable specifies the part of the stress tensor from 
which external work can be extracted. Another of the hydrodynamic variables is the velocity 4-vector (various 
choices of u are discussed in d and (H). Finally one has the dissipative part of the stress tensor. When bulk viscosity 
is neglected, this is a traceless symmetric tensor, tt^''. Since this part expresses shear, it must vanish on contraction 
with u^. 

Since we will discuss longitudinal flow, it is convenient to go from Galilean coordinates (i, x, y, z) to the curvilinear 
coordinates (r, rj, r, 0) where t = — z^, ^ = tanh^^(z/t), r = \/ -V and = tan^^(y/a;). The metric becomes 
gp_y = diag(l, — T^, — 1, — r^). The only non-vanishing Christoffel symbols are FJJ^ = r, F^^ = — r, F^^ — F^^ = 1/r 

and F^^ = F^^ = 1/r We will write covariant derivatives as and partial derivatives as d^. The action of on 

a scalar fleld is the same as the action of 9^^. In terms of the proper time T — Vt^ ~ r"^ , one deflnes the components 
of the velocity 4-vector — dx^ /dT. One can show that u^^u^ — 1. 

Longitudinal flow is an approximation applicable to ultra-relativistic heavy-ion collisions when the hydrodynamic 
variables at any point of spacetime depend on r and rj but not on r and <f>. Clearly such an approximation is valid far 
from the edges of the fluid volume, and at times r ^ R/cg, where i? is a typical transverse size. We will parametrize 
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the velocity vector by a quantity y in the form — 

u'^ = (coshy, — sinhy, 0, 0). (1) 
r 

We note that the rapidity is y + 77. Scahng flow corresponds to y = in our notation, as we discuss later. A 
fluid element with y = 0, in our notation, moves along a world line of fixed 77, corresponding to a constant velocity 
V = tanhy^. The divergence of u is 

6 = dnu'' = Ur sinhy + (y„ + 1)- cosh y. (2) 

T 

This defines a macroscopic time scale for a hydrodynamic flow. Here, and later, we use the notation Jt etc., to denote 
the derivative of a scalar / with respect to the variable t etc. We also define the material derivative D = u'^dfj., and 
through it the spacelike vector Du'^ and its norm = —Du'^'^Du^. A straightforward computation shows 

S = yr cosh y + (yrj + I) - sinh y. (3) 
r 

Then it is easy to write down the unit space- like vector = {Du'^)/S, 

= (sinh y, - cosh y , 0, 0) . (4) 

T 

One also defines another directional derivative operator D — v^^d^. In the local rest frame one finds that D is the 
derivative with respect to time and D is the longitudinal spatial derivative. 

Using the methods outlined in the appendices, one finds the hydrodynamic equations — 

De + BQe = GTTy, 

clDe + BSe = Dnv + Snv, 
4 

t^Dttv+ttv = -^VvQ, (5) 

where B ~ 1 + c^, jjy is the coefficient of shear viscosity and r^r is the relaxation time associated with the shear stress. 
Try. Recall that such a relaxation time is necessary to construct causal hydrodynamics jisl . [l^ . For longitudinal 
flow we have reduced the tensor equations of hydrodynamics to three coupled scalar equations for the three scalar 
hydrodynamic quantities e, y and Try. 

ELNS hydrodynamics is recovered when = 0, so that the last of eqs. ([5]) reduces to Try = 477^0/3. One expects 
that for fluids which evolve inertially, i.e., in the absence of external forces acting during the evolution, the solutions of 
eqs. ([5]) should approach the solutions of ELNS hydrodynamics at times r 3> tv. This conclusion may clearly change 
when a fluid is acted upon by external forces at all times. In such cases, of course, driving terms have to be added to 
the equations. 



A. Material properties at vanishing chemical potential 

For an ideal fluid, the equation of state can be cast into the form p = c^e, where p, e and Cg, could all be functions 
of the temperature, T. Straightforward dimensional analysis shows that e — bT'^, where b is dimensionless. In general 
there are various intrinsic mass scales, fj,i, in the fluid, and b could have an implicit dependence on T through the 
functional dependence, b{T / fjLi,T / ^2, ■ ■ •)• If '^s = 1/3 at all temperatures, then the trace of the stress tensor vanishes 
identically. This implies a special symmetry called scale symmetry, or conformal symmetry [23|. One aspect of 
conformal symmetry is that uniform scaling of external scales such as T by a constant leaves material properties 
unchanged. Clearly, then b cannot depend on T, and must be constant. 

For a non- ideal fluid, the stress tensor has an additional viscous part. When the trace of the full stress tensor 
vanishes, then the bulk viscosity vanishes identically. Conformal symmetry implies that the physics of such fluids 
can be expressed in terms of dimensionless combinations of material properties which are temperature independent. 
One such combination which has been used in the literature is rj^/s, where ry^ is the coefficient of shear viscosity 
and s — {e + p)/T is the entropy density. In passing we note that for a conformal fluid, s = je^^^, where 7 is a 
dimensionless constant. 

Causal viscous hydrodynamics requires another material property of the fluid, the relaxation time for the shear 
part of the viscous stress tensor, Tj^. Dimensionally, = a/T, where a is dimensionless, and becomes constant when 
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the fluid has conformal symmetry. This dimensionless number is proportional to the quantity called liquidity [28| 
which, in non-relativistic fluids, measures the mean-free path (proportional to r^) in units of the interparticle spacing 
(proportional to 1/T). In a gas, this number is very large, in liquids, small. It follows from the expressions 

for s, T^r and the equation of state, that s = /CeT^r, where the dimensionless constant /C = 4/3a. This implies that the 
dimensionless material property eTT^/rjy = s/JCr]v We call this combination x and discuss it extensively in the next 
subsection. 

In reality, the fluids that we are interested in are not conformal . At temperatures below Tc the fluid of hadron 
resonances has a plethora of mass scales, which breaks conformal symmetry, and manifests itself in deviations of cj. 
from the value 1/3. It is not a big stretch of the imagination to expect that bulk viscosity will be non- vanishing in this 
fluid. A fluid of quarks and gluons also breaks conformal symmetry through the conformal anomaly, which results 
in the running of the strong coupling and the appearance of the QCD scale Aq^d, and explicitly through the quark 
masses. One question of interest is how important are these departures from conformal symmetry. 

One might expect that at very large T, when the QCD coupling is close to zero, and all the quark masses are 
much less than T, one might have conformal symmetry to a good approximation, by virtue of the fluid being well- 
approximated by a massless ideal gas. In fact, this is the limit in which the Boltzmann fluid approximation is seen to 
hold in weak coupling theory, with x = 3/2c^ = 9/2. In this limit, one has, additionally, bulk viscosity much smaller 
than shear viscosity [30|, ~ 1/3, and the energy density close to an ideal gas value. 

Lattice computations show, surprisingly, that approximate conformal symmetry is obtained also at substantially 
smaller T/T^ «2-3, where the pressure deviates significantly from its ideal gas value [31j. Toy models of QCD with 
substantially enhanced {N — 4 super-) symmetries, which give up the running of the coupling and asymptotic freedom, 
have been used to model this observation. They are bound to fail in the vicinity of where the conformal measure 
[3H is large, and bulk viscosity cannot be neglected [l^l- One prediction from these toy models, using the AdS/CFT 
conjecture, is that r]v/s = l/47r, yielding x = 47r//C. A recent computation in an appropriate A'^ = 4 SYM theory has 
found a = (1 - ln2)/67r [H, which then yields x = (1 - ln2)/2 w 0.15. 

In this paper we shall examine three models of viscous fluids. The flrst, which we call a simple fluid, is one in 
which the material properties c^, rjv and are constant. Lattice results show that Cs is almost constant over a range 
of T 31]. However, preliminary lattice computations of transport coefficients are almost consistent with the power 
counting in T over the same range {34|. As a result, the main motivation to study this model of a simple fluid is 
not its direct application to heavy- ion physics, but the fact that it allows explicit computation of the hydrodynamics, 
and contains qualitatively all the phenomena that we find with other models of viscous fluids, as we show in a later 
section. 

The more restricted models of materials that we use have the property that x is constant. A Boltzmann fluid 
has been examined in the literature [l3|, and is defined by the specific value x = 9/2- We devote one section to 
detailed hydrodynamics of the Boltzmann fluid. In addition, in a subsequent section, we examine the whole class of 
conformal fluids with x — 37ra, for various a. Note that the hydrodynamics of the conformal fluid with a — 3/27r 
(i.e., — 3/2ttT) is exactly equivalent to that for a Boltzmann fluid. As a result, it does not seem possible to use 
hydrodynamics alone to distinguish a conformal fluid from a Boltzmann fluid. We discuss this in greater detail later. 



B. Laws of physical similarity 

The equations of ideal hydrodynamics are obtained by setting ttv = rj^ = in the eqs. Then the equations for 
the remaining hydrodynamic variables, e and y, are — 

De + BOe = 0, c^De + BSe = 0. (6) 

Consider the symmetries of these equations. The solutions are unchanged by the independent scalings e Ae and 
r —^ Ct. Introduce the variables e = ln(e/eo) and 9 = ln(r/To), where the arbitrary scales eo and tq can be chosen to 
be the initial conditions. This removes the freedom of scaling, so that the solutions of the equations can be written in 
the form e{9,ri) and y{9,ri). Thus, the scale symmetries of the ideal hydrodynamic equations connect solutions with 
different initial conditions. 

Non-ideal hydrodynamics breaks both these symmetries by the introduction of the time scale and the scale of 
energy density vj — rjv /t,^. In other words, the scaling e ^ Ae (simultaneously Try Attv) and t C,t are not 
symmetries unless one simultaneously scales — *■ and rjv — > A^r^v-- Thus, the scalings relate flows of fluids 
with different material properties. This is the relativistic analogue of scaling laws called "physical similarity" 3] that 
one finds in non-relativistic fluids. Such similarities are the basis of scaling invariants, also known as dimensionless 
variables, such as the Reynolds number. Re, which are used to relate flows of different fluids. 
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The analysis here gives three scahng invariants — 

X=— , ^ = — , and S = — . (7) 

WW e 

The dimensionless ratios lead to physical similarities between flows. We can relate these variables with quantities 
familiar from Navier-Stokes hydrodynamics by examining what they become in the appropriate limit. 

As pointed out earher, the ELNS hmit of eqs. (O is obtained when = 0. In that case, Try = 4ryv6/3. The 
quantity is the inverse of a characteristic scale for the flow, Tj,. In the non-relativistic, i.e., the Navier-Stokes, limit, 
Tc is a characteristic time scale. In this limit we can define a characteristic length scale for the flow by the relation 
i-ic = TcV^ where v is the flow velocity. Then, one finds 

477V ^4c|/^y^^4c|M2 
3eTc 3 \Cs) eLcV 3 Re ' 

where M — v/cs is the Mach number of the flow and Re — eLcv/rjv is the Reynolds number. The first of the 
expressions on the right comes from taking the ELNS limit, whereas the last expression involves taking, additionally, 
the non-relativistic limit. In Navier-Stokes hydrodynamics the ratio of inertial and viscous forces is Re In this 
case one may therefore suspect that S cx 1/Re. The exact relation above bears this out, with corrections needed to 
translate between the fully relativistic and non-rclativistic formulas. Similarly, one finds that 

where a mean-free path, A = ^.CsTtt, S, is some numerical constant, and K = X/Lc is the Knudsen number. 
The third variable 
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is interesting, since the combination KKe/M ~ 1 in Navier-Stokes theory. For Boltzmann and conformal fiuids this 
combination is constant. The Navier-Stokes relation is obtained for these fiuids when ^ is chosen appropriately. For 
the Boltzmann fiuid, this happens when ^ ~ 2/3. Eqs. ((8l [9l fTO|) provide the connection between S, (p and x and Re, 
M and K in the appropriate limit. 

In ideal hydrodynamics thermalization and freezeout are notions which are imposed from the outside. In non-ideal 
hydrodynamics, however, some understanding of these phenomena could be possible (ssj . In the Navier-Stokes theory, 
for example, K ~ Af/Re must be smaller than unity in order for the solutions to describe valid flows. Since a solution 
of the Navier-Stokes equation allows us to compute both M and Re, one can use the solution to compute K and 
determine its own validity. A solution of the equations of IS hydrodynamics gives Tc — 1/8. When Tc is larger than 
T,r, i-G., when (p < 4/3, the solution corresponds to a physical flow. In the case of scaling flow (discussed below), this 
gives an initial time, tq, at which the solutions begin to describe physical fluid flows. Thus we have a self-consistent 
description of thermalization. For scaling flow, ip decreases with time. Hence, after thermalization, IS hydrodynamics 
is always applicable; i.e., we lack a description of freezeout. That phenomenon requires us to examine radial flow. 
As argued before, radial flow becomes important at a time tt — R/cs, where R is the transverse size. Thus, one- 
dimensional IS hydrodynamics in the scaling approximation is expected to be valid in the range tq < t < tt- A 
description of freezeout at late times has to be sought in the full 3-dimensional hydrodynamics. 



C. Scaling solutions and fluctuations 



Solutions with y = are called scaling solutions or boost-invariant solutions. It was argued by Bjorken that 
asymptotic freedom implies that, at sufficiently high energies, hadron multiplicities must become invariant under 
longitudinal boosts. Assuming further that these multiplicities have their origin in the hydrodynamic distribution of 
the entropy density, he argued that the relevant fiows in high energy heavy- ion collisions must be boost invariant [5]. 
While the phenomenological relevance of this argument may be questioned, boost invariance is simply analyzed. One 
uses ?/ = 0, as a consequence of which S = and O = 1/t. Substituting these into eqs. ([5]), one finds that 

~ 2 ~ T 4 
tDc = ttv — Be, Dttv = c^De, tDtty = Try -I- -zu, (11) 

Ttt 3 

where D = dj- and D = (1/t)(9^ (see eq. IA4p . Now, Fourier transforming in 77 decouples the Fourier modes, labeled 
by k. The second equation implies that the identity 7ry(r, k) = c^e(T, k) must hold for all fc 7^ 0. However the other 
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two equations cannot be manipulated to give D(7ry/e) = 0. Consequently, only the k = mode is allowed to be 
non-zero, i.e., the solutions to these three equations must have both e and ny independent of rj. This demonstrates the 
well-known equivalence of the conditions of boost invariance and y = 0. Of course, the second of the three equations 
above becomes redundant and the problem can be treated with the remaining two equations. 

In the next few sections we will investigate the scaling solutions for fluids with various different constitutive equa- 
tions, and analyze their stability 0, [1^. Some numerical studies of the correlations of fluctuations were reported in 
[2^. We note that the IS hydrodynamic equations can be written in the form xg = f(0,x, x^), where x is the vector 
of three hydrodynamic variables and 9 is the "time" variable. Now, setting to zero the component corresponding to y 
in X results in the boost-invariant equations, whose solutions we represent by x°. Represent the fluctuations around 
this solution by Ax-*^, where A ^ 1 is a dimensionless book-keeping parameter. The equations for these fluctuations 
can be written to linear order in A, 

xi = Mo(0,x")xi + Afi(0,xO)x;, (12) 

where Mq and Mi are the Jacobian matrices of the derivatives of f with respect to x and x^ respectively. For the 
stability analysis one asks whether a given x^ increases or decreases with time. Fourier transforming in t] decouples the 
derivatives with respect to the variables and gives independent linear evolution equations for each mode — Xg(0, k) — 
Mx^{9, k), where M — Mo +ikAIi. The question of stability then reduces to examining M and checking whether the 
solutions for x^ decrease faster than the scaling solution or not. We demonstrate the method with the ideal fluid in 
the next section. 



III. IDEAL FLUID: THE BJORKEN SOLUTION AND SOUND WAVES 




k k 



FIG. 1: (Color online) On the left is the dispersion relation for sound waves in the Bjorken solution for — 1/3. The real 
part of u is shown in blue as a function of k and the damping exponent, i.e., the imaginary part, is in red (the straight line in 
black is the line a; = Csk). There are no propagating waves for k < l/y/S. On the right is the power spectrum resulting from 
the evolution of Pe{0 = 0; k) = 1 a,t 9 — 1 (upper curve) and 3/2 (lower curve). 

The scaling solution in ideal hydrodynamics is obtained by setting rjv = ttv = in eqs. pTjl . Introducing the 
variables e = ln(e/eo) and = ln(r/ro), the equations become — 

de , de ^ , ^ 

This gives the Bjorken solution 



ei9,v)^-B9, i.e., e(r, r?) = eo (^) . (14) 



The entropy density is defined as s = (e + p)/T. Entropy is conserved in an ideal fluid 0]. All this is well known. 

The linearized equations for sound waves are extracted by introducing a dimensionless parameter A <C 1, which sets 
the scale of quantities involved in the propagation of sound relative to the boost-invariant background, as discussed 
earlier. Then one proceeds by setting y{0,ri) — Ay^{6,ri) and writing e{9,rf) = e^{9) + Ae^{9,r]). Inserting these into 
the ideal fluid equations ([6]), one finds to order A", the equations p^ . At order A one finds — 

l^-Byl^, and yl = - (l - ^\ e^ + {B - 2)y\ (15) 



e 
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A Fourier transformation in 77 reduces the partial differential equations above to uncoupled evolution equations for 
each Fourier mode. Since these linear first-order evolution equations are autonomous, i.e., they do not involve the 
variable 9 anywhere explicitly, the solutions are in the form of waves exp[i(w0 + fcry)]. 

Before proceeding, a point about wave solutions may need comment. Plane waves in Galilean coordinates have the 
form exp[i(a;i + k • x)], unlike the one above. However, the functional form of a wave solution is dependent on the 
geometry of the situation; for example, in spherical geometry a wave solution has the form cxp[i(Ljt + A;|x|)]/|x|. Wave 
solutions in boost-invariant geometry have the specific form 



A real value for lu corresponds to an oscillatory solution. In the form on the right it is manifest that the real parts of 
Lo and k are not inverses of typical time and length scales for oscillation. If On and are, respectively, the n-th and 
n — 1st times that the phase of the wave repeats at a given spacetime rapidity, then clearly On — On-i = 2tt/uj, so that 
Tn = Tra-i exp(27r/ti;). In other words, the "period" of oscillation increases geometrically in the number of periods. 
The n-th time that the phase recurs after the initial time tq is given by r„ = tq exp(27rn/ci;). In the same way, at fixed 
t, the "wavelength" increases with z. Both these scalings are direct consequences of boost-invariant expansion — the 
longitudinal expansion red shifts sound waves. The analogy with the red shifting of light in an expanding universe is 
clear; both follow from the fact that the spatial components of the metric depend on the time. If lu is complex, then 
the real part gives oscillations exactly as described above. The imaginary part of uj, i.e., — Reiw, gives rise to power 
law behaviour in r. The scaling solution is stable when Reiw < 0, so that fluctuations are damped. Hence we will 
give the name damping exponent to — Reiw. 

Substituting the form of the wave solution in eq. (|16p into eq. psp . one obtains the dispersion relations 

*'^ = -^(l-a±iy^4c2(fc2_fc2) ^here = (17) 

For k < ko the modes are overdamped; the two damping exponents are equally spaced around (1 — c1)/2. Only for 
k > ko, are there propagating modes. These are damped due to the expansion of the scaling solution. The slowest 
decreasing part corresponds to the positive sign above. As a result, this is the part that is visible to experiments. 
This dispersion relation is shown in Figure [T] A general solution for can be written in the form 

e\0,k) ^c+e"^+'' + c-e"^-^, (18) 

where c± depend on the initial conditions, and uj± arc the solutions in eq. ()17|) with the corresponding signs. Similar 
solutions can be written for y^. Since the real parts of iu) are non-positive, the fluctuations do not grow, and the 
scaling solution is stable [2^ • 

The physics of sound can be captured in the evolution of the power spectrum of fluctuations of the energy density, 

P,{T;k)^\e\T;k)\^ where e\T;k) = f ^e~''''Je\T,7j). (19) 

J V 27r 

Since e^{T,ri) = e^{T,ri)e^{T), for the ideal fluid one may write asymptotically, when the component in c_ can be 
neglected, 

Pe(T;fc) =P.(To;fc) (^)'"', where Wrf(fc) = 3 + ^(1- c^) - Cs^k^ - fc2H(fco - k), (20) 

ko is given in eq. (jl7p . and H denotes the unit step function. At not so late times, the interference between the 
frequencies and uj_ (when k > ko) gives rise to beats. The shape of the power spectrum resulting from an initially 
flat power spectrum of fluctuations through exact solution of eqs. (jisp with initial conditions c+ = c_ is shown in 
Figure [TJ At short times it is dominated by beats. The expression in eq. is asymptotic. As shown in Figure [1] it 
is not recovered for — 3/2. It is clear that if the initial conditions contain fluctuations around the boost-invariant 
values, then after sufhcient evolution, these fluctuations are damped. The longer the duration of hydrodynamic 
evolution, the more damped the fluctuations. A detailed analysis of the growth of transients at short times is given 
in Appendix C. 

If the initial conditions, i.e., the power spectrum at time tq, for an ideal fluid were known, then an event-by-event 
measurement of the power spectrum of the acoustic energy density would be able extract the value of and thereby 
give a measurement of the equation of state. Conversely, if the equation of state were known, then the same data 
could be used to extract, event by event, the initial conditions. 
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The power spectrum studied here is closely related to the correlation function of fluctuations studied in ^2^. In 
Fourier space the correlation function corresponds to studying the joint distribution of fluctuations at different fc, 
whereas the power spectrum gives the variance in the fluctuations at a single k. Sonic peaks are visible in both the 
quantities. The analysis of the power spectrum in terms of interference, and the consequent clear relation with c^, is 
a little harder to establish for the correlation function. 



IV. A SIMPLE FLUID 



The simple non- ideal fluid model defined in Section II. A has constant values of all constitutive parameters — Cg, 
and rjv Using the variables x a-nd (p introduced in eq. ([7|) and z9 = t/ttt, we make the decomposition — 

xi^,v)^X°{^) + ^xH^,v), ^(^?,77) = ^°W + A'/''('?,^), y{^,v) = ^y\^,v)- (21) 

These expansions have to be substituted into the eqs. ^ and the material properties of the simple fluid used to 
extract equations for the boost-invariant solution x° and ip'^, and the fluctuations x^j 3i,nd ip^. In the next two 
subsections we examine these two problems. 



A. The scaling solution 



The equations satisfied by the scaling solution are- 



and 



One solves the second equation and then inserts the solution into the first. It is easy to check that 



6+-Ei(^) 



& = --Ei(^o)+e'^V"(i^o), 



(22) 



(23) 



where Ei(x) denotes the exponential integral [3a |. and i9o — tq/t.^. Wc will use i^o = 1 in numerical work. The 
asymptotic expansion of the exponential integral, 



Ei(a;) 



e=^ / 1 
- 1 + - 

X \ X 



can be used to write down the asymptotic expression — 



(24) 



(25) 



Note that the solution of the homogeneous equation decays much faster than the particular integral. Hence, the 
long-time behaviour of Try is nearly independent of the initial conditions on this quantity. 




FIG. 2: (Color online) Comparison of ELNS hydrodynamics (lines in black) and IS hydrodynamics (in red). The evolution of 
normalized energy density, x", normalized shear stress, and S = •pP jx^ , are shown as a function of the normalized time, 
t/tt,. The initial condition ^^{V) ~ 4/3 in both cases. 
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Inserting the solution in eq. (|23|) into the first of eq. ((22)) gives a solution in terms of the Meijer G-functions (36l.[37t. 
Instead, by inserting the asymptotic expansion in eq. (|25p into the equation, we can find a simplified solution — 

xV)-^-^V^ + 3^ + --- (26) 

where a is determined by the initial condition x*^(i9o) and the remaining terms decay as higher integer powers of l/i?. 
The incomplete Gamma function is defined to have the values T{B, 0) = r(_B) and r{B, oo) = 0. The solution of the 
homogeneous equation gives the Bjorken solution. However, this is not the slowest falling part; that is given by the 
l/'d term, which arises from the inhomogeneous part of the equation, and has no dependence on initial conditions. In 
the very late time limit, one finds (p'~' ~ 4/3?? and x° ~ 4/3c^??, so that the Bjorken solution is never recovered. The 
expansion is slowed down with respect to the ideal fluid because of dissipation — in expanding against the vacuum, 
an ideal fluid pumps all its energy into the expansion, whereas a non-ideal fluid dissipates energy, thus slowing the 
expansion. Similar results were also obtained in ^ for the case of ELNS hydrodynamics. 

1. Comparison of ELNS and IS hydrodynamics 

The differences between ELNS and IS hydrodynamics are illuminating. The ELNS limit of IS hydrodynamics is 
obtained formally by setting = 0, or, equivalently, by dropping the term Dttv in eqs. ([5]). Now, a comparison of the 
two theories makes sense only after appropriate normalization of both time and energy. So in ELNS hydrodynamics 
one must introduce artificially a unit of time which is numerically equal to the value of t^^ used in IS hydrodynamics. 

Then, rewriting the ELNS equations, one obtains scaled quantities and (p^ which are directly comparable to the 
same quantities in IS. The equations of ELNS are then obtained from eqs. ([^ by dropping the term dip^/d'd, so that 
ip'^ = 4/319. Inserting this into the equation for one has 

A comparison of the resulting solutions is shown in Figure [21 The initial conditions for ip*^ are chosen to be the 
same; in IS hydrodynamics it corresponds to the smallest time at which the equations are valid, i.e., tq = Tjt. As in 
[ll|,[lj], we find clear differences between ELNS and IS evolution. In earlier works the difference in the evolution of 
the energy density was investigated for equal initial conditions. We have chosen the initial conditions for x° to be 
different, but tuned so that the ELNS and IS solutions approach each other at large times. Since initial conditions are 
not detectable in heavy- ion collisions, it is of interest to see that large-time behaviour cannot, by itself, distinguish 
between ELNS and IS hydrodynamics. Moreover, the same final energy density can lead to different estimates of the 
initial energy density in the two kinds of dynamics. 

2. Entropy production 

Assuming that s = 76'^^'*, where 7 is a dimensionless quantity which is almost temperature independent, a = (x'')'^^^ 
is a dimensionless quantity proportional to the entropy density. The first of eqs. (j22p can be easily manipulated into 
the form — 

'dtd " 4i9ai/3- ^ > 

Using the solutions for (pP and x" in eqs. ([25l [26]) . one finds that at late times the right hand side is negative. This 
is in accord with the exact solution for x° shown in Figure O Depending on the initial conditions, a may increase 
initially. However, at sufficiently large time it must decrease. If there is initial growth in a, then the turnover comes 
when the right hand side of the above equation passes through zero, i.e., at the value of d when x° — 'i^p^ /A. 

Note, however, that the element of 3- volume contains a scale factor r from the metric. This implies that the total 
entropy scales as S = cri?. It is a straightforward exercise to rewrite eq. (^5)) to obtain 



3^° / 19 ^ ^'^ 



d-d 4 V S 



(29) 



The right hand side is manifestly positive definite, indicating that the total entropy increases with time. Using the 
asymptotic expansion of tp^ in eq. (j25p . we find that E grows asymptotically as i?^/"*. One has the same power law 
growth of S in ELNS dynamics. 
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B. Sound waves 



Using the decomposition of cq. (j2ip . we examine small fluctuations around the scaling solution. At large 1} we may 
use the asymptotic solutions x° = 4/30^19 and (p'^ = 4/3^? in the fluctuation equations. We also transform to the 
variable = Ini? and use = gexpd. After Fourier transforming in 77, the equations for fluctuations take the form — 



X 



M 



Ml = 



/ 



3c: 
4 

V 



where M = ikMi + Mq 



4 

3c? 



Mn = 




(30) 



In contrast to the equations for fluctuations in an ideal fluid, i.e., eqs. (jl5p . these equations are not autonomous. As a 
result, they cannot be solved by Fourier expansion in 6. A numerical solution is always possible, and we can examine 
the limits of large and small k analytically. At every r there is an upper cutoff on k imposed by the requirement of 
the applicability of hydrodynamics. This cutoff increases with r. 




0.2 0.1 0.6 0.8 



1 1.2 1.4 




FIG. 3: (Color online) Solutions of the linearized equations for fluctuations around the scaling solution in a simple fluid. The 
normalized energy, y shear, (p^, and are shown for k — 1/2 (red), 1 (green), 2 (blue), 4 (purple) and 8 (orange). The flrst 
three values of k are overdamped, but the last two show oscillatory behaviour. The frequency of oscillation increases with k. 




FIG. 4: (Color online) The panel on the left shows Inx^ {k,r) as a function of 6 = ln(r/r^) for the simple fluid with — 1/3 
for fc = 1/2 (red), 1 (green), 2 (blue), 4 (purple), 8 (light blue) and 16 (orange). The panel on the right shows the damping 
exponent, i.e., the imaginary part of u, for — 1/5 and 1/3, plotted as a function of k. 

In the limit fc — > 0, one may set M — Mq. Then, since A/g is diagonal, one can read off the solutions easily, s-nd 
g (and hence y^) are overdamped, whereas ip^ decays exponentially in r. The exact solutions are — 



y\k^O,T) = yi(fc^O,ro)(^), 
<^i(fc^O,r) = <^i(fc^O,To)exp(^^i— . 



F B 



r B 



To 



(31) 
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In the regime fc ^ 1, a first approximation would be to neglect Mq. As a result, one expects uj (x k, where iui is an 
eigenvalue of M. Since w is very large, exp0 changes little over many oscillations. Consequently one could treat this 
factor as constant whenever it appears inside M. Within this approximation the equations above can be treated as 
autonomous and therefore generically describe oscillations. 

1 




2 4 6 8 10 12 14 
k 



FIG. 5: (Color online) The power spectrum obtained from an initially flat power spectrum of fluctuations at ^ = 1 (red), 1.1 
(blue) and 1.5 (black). Here — 1/3 and tq has been set equal to r^. Initial conditions are X^(0) ~ 3/^(0) ~ V^(0) = 1. 

Since the eigenvalues of Mi are real, within the approximation where one neglects Mq, the solutions are not damped. 
To go beyond this and obtain the damping exponent, one sets up a perturbation theory in powers of l/k by writing 
iuj — ikXi + Aq + 0{l/k). The simplest way to proceed is to substitute this form into the characteristic equation for 
M. One finds then that this equation has a leading term of order k^ (which solves the eigenvalue problem for Mi) 
and a first perturbation term of order fc^. The latter shows growing solutions in the unphysical region exp0 ^ 1, but 
damping with iiv = —2 + 0{l/k) when exp9 » 1. Note that this damping exponent is independent of c^. 

The numerical solutions to eqs. ([50]) are shown for a range of k in Figure [31 There is clear evidence of overdamped 
solutions for fc < 2 and damped oscillatory solutions for fc > 2. In Figured] we show that x^, which is proportional to 
the energy density of fluctuations, is damped fairly rapidly. However, for < 1 there are clear signs of transients; a 
detailed discussion of which is given in Appendix C. Numerically we see that is damped as a power of r at large 
T, making it easy to extract the damping exponent numerically. As shown in the figure, and in agreement with our 
analysis above, at small k this goes to 1 + and climbs to the neighbourhood of 2 at large fc. 

The power spectrum of fluctuations starting from an initially flat spectrum, Pc(0;fc) = 1, is easily amenable to 
computation, and is shown in Figure [5] The transient growth phenomenon gives rise to several peaks for < 1. By 
^ ~ 1.5 the effect of the damping exponent is clearly visible. A complete analysis is given in Appendix C. 



1. Sound in ELNS hydrodynamics 

Sound propagation in ELNS hydrodynamics is quite different. The equations for sound can be obtained from eqs. 
([SO)) , by simply dropping the term dtp^ /d9, i.e., by setting iy9^ — Aikg/3, and treating as an arbitrary scale of time, 
set equal to the value it would have in IS hydrodynamics. When this relation is inserted into the equation for gg, a 
term in g,;^ is generated. This is the diffusive term we expect in ELNS hydrodynamics, and arises directly from the 
neglect of the relaxation time in shear pressure. 

As before, we transform to variables = lnz9 and write = gexpS. The equations for sound in ELNS hydrody- 
namics are 

dg(^^^=M(^^\ where M = - [ + | . (32) 

\9j \9j' \^|icffc 2 + clP J 

In the limit of small k, i.e., fc ^ 1, it is clear that oc t^^ and oc l/r. This is similar to the results in eq. (|3ip . 
In this limit the solutions of ELNS and IS hydrodynamics are not qualitatively different. This is expected since the 
characteristic time scale of the wave is much larger than r^r. 
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FIG. 6: (Color online) The power spectrum obtained in ELNS hydrodynamics for a simple fluid from an initially flat power 
spectrum of fluctuations at 6 — 2 (red) and 2.5 (blue). Here = 1/3 and to has been set equal to Tn. 



In this case, an exact solution is possible for all k, since the equations are autonomous. The eigenvalues of M are 

where the combinations 4(1 — c^)x^/3Cg — i/c(7 and Aikx^ /Scl + g decay respectively with these damping exponents. Not 
only is diffusively damped at large k, but there are no propagating modes at all. This behaviour is characteristic 
of parabolic equations. 

In the absence of propagating modes there are no beats. Transient growth can occur, but there is only a single 
peak in the power spectrum of the fluctuations of energy. The transient analysis is given in Appendix C. The power 
spectrum resulting from an exact numerical solution, starting from Pe(0; fc) = 1, is shown in FigurelHl Comparing this 
with Figure[ni we see that there is a clear difference between diffusive damping of fluctuations in ELNS hydrodynamics 
and sound in IS dynamics. 



V. A BOLTZMANN FLUID 

A Boltzmann fluid is defined by the constitutive relation x — ^'''■rr/Vv — 3/2Cj — 9/2, where, as discussed in Section 
II. A, the three quantities e, ?7v and r^r, all depend on the temperature. We rewrite the hydrodynamic equations in terms 
of the variables S, y, u = Tt / {TqTo) , and = ln(r/To), where the initial conditions are applied at the time tq, i.e., at 
6 = 0. We analyze the scaling solution and its stability by the usual technique of writing u{t, rf) — u^{t) + /S.u^{t, j]), 
S(r, 77) = S°(t) + AS"'^(t, 77) and y{T,r]) = Ay^(T, 77). Substituting these into eqs. ([5|), using the material properties, 
and separating out the equations to orders A" and A, we obtain the equations which lead to the scaling solution from 
the former, and the equations for fluctuations from the latter. 



A. The scaling solution 



The equations for the scaling flow become 



V, 



where 



V 



u [S + 3 - c2] /4 
- S2 + S [(1 + c2) - Bu] 



(34) 



A = 4/3%, B = To/r^(0) and are non-negative. In this subsection we lighten the notation by writing u for and 
S for S°. Three numbers are needed to flx the initial condition in the original formulation of the problem (eq. [5|), i-e., 
the initial time tq and the values of e(To) and 7ry(ro). The initial condition on u is, by definition, u{9 = 0) = 1. The 
two free parameters in the initial conditions are transformed into the value of S(0) and the value of the parameter 
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B which appears in the equation. Note that S < 1 is disfavoured. A dimensional quantity is needed to complete 
the specification of the initial conditions, and we can choose this to be Tq = qB/tq. Next, choosing = 1/3 (and 
hence x = 9/2), one has A = 8/27. Since the equations are autonomous, one can analyze them using the phase plane 
method [H]. 



1. Phase plane structure 




FIG. 7: (Color online) The panel on the left shows structural elements of the phase plane: the nuUclines, i.e., the lines along 
which one of the derivatives vanishes (that of it vanishes on 770 and ry, , that of S vanishes on (t±), and the fixed points, at which 
both derivatives vanish. The panel on the right shows the vector field of derivatives and the fixed points. 



The idea of the phase plane method is that the right hand sides of eqs. define the local direction of derivatives, 
which is uniquely given by coordinates (u, S) on the phase plane. Knowing this vector field, the solutions of the 
equations are simply integral curves obtained by following the vector field starting from any initial condition. A 
special role is played by the zeroes of this vector field, i.e., stable solutions or fixed points of the differential equations, 
and by nuUclines, which are the lines along which one or the other derivative vanishes. 

The nuUclines of u are the lines 770, which corresponds to w = 0, and 77*, which is S = — 3. The nuUcline 770 
happens to be an integral curve, since the vector field is everywhere tangent to this line. The nuUclines of S are the 
hyperbolae + BS[u — {1 + c'^)/B] ~ A. This has two branches, cr_ has S < and is asymptotic to the u-axis at 
u = —00, tT+ has S > and is asymptotic to the M-axis at u = 00. 

These two sets of nuUclines cross at three fixed points — P+ = (w+,S+) is the intersection of rjo and (T-|_, P_ = 
(u_, S_), which is the intersection of 770 and (t_ and = (u^, S*) which lies on 7/* and cr_, and 



u± = 0, 



1 

B 



A 



l + 4± yJ^A + (1 + c2)2 



S, = 



(35) 



These features are shown in Figure [7] 

In a small interval around any fixed point (li^, S/) one can linearize the equations to get 



M(u/,S/)[g:g^J, where M(7i, S) - ^ ^ + 



u/4 
-2S + (l + c?) 



Bu 



(36) 



At P± one of the off-diagonal components vanishes as a result of which one can write down the eigenvalues by 
inspection. P-|_ is a hyperbolic fixed point. P_ is a repulsive fixed point for a Boltzmann fluid, but changes into a 
hyperbolic point for large enough A. One also finds that P* is a hyperbolic point for a Boltzmann fluid but changes 
into a repulsive fixed point for sufficiently large A. The reason for this is not hard to understand. At ^ = 4(3 — c^) 
the points P_ and P* are coincident and exchange character, leading to a saddlc-nodc bifurcation at this value of A. 
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FIG. 8: (Color online) The solutions corresponding to the unstable manifold for the normalized energy, e/eo = u(^)*e~^*, the 
quantity S, and the normalized shear stress, Su(0)''e~*^, plotted against r jT-^. In the figure for the normalized energy, the 
curve in yellow shows the Bjorken solution normalized to match the exact numerical solution at large r/r-n. 



2. The unstable manifold of P+ : a physically relevant solution 



We concentrate on the physical flows, i.e., those originating on the half-line u = 1 and having S > 0. Since the 
vector V points into the part of the plane bounded by m = 1 and S = 0, one sees that physical flows remain in this 
portion of the phase plane. Also, since there are no fixed points in this part of the plane, one finds that flows cover 
the whole quadrant. Using Bendixson's theorems [39j . one can prove the intuitively obvious fact that there are no 
periodic or quasi-periodic solutions starting from these initial conditions. 

The main organizing element behind the physical flows is a special solution, U, i.e., a curve in the phase plane, 
called the unstable manifold of P+. This is the solution with initial conditions in the infinitesimal neighbourhood 
of P+, the explicit numerical solution for which is exhibited in Figure [H] Since P+ is a hyperbolic fixed point, all 
solutions starting from initial conditions above U are attracted to it from above, and those starting below it are 
attracted towards it from below. Clearly, then, an important element of the analysis is Su, the point at which U 
intersects the line of physical initial conditions u = 1. 

Two varieties of stability in the solutions should be noted. At late times the flows are stable against changes in 
initial conditions, since they are always attracted to U. One useful consequence is that fairly crude PDE solvers suffice 
to integrate eqs. (IM]) . There is a deeper level of stability, called structural stability, arising from the fact that P+ does 
not collide with any of the other fixed points on changing the parameters cj., A, and B. This means that for all values 
of these parameters, the nature of physical flows is determined by the corresponding solution for U. The usefulness 
of structural stability is that extraction of parameters from data becomes particularly simple. 

The first step to solving for the flows is to determine U. In the neighbourhood of P+ it coincides with the 
eigenvector corresponding to the positive eigenvalue of M(m+, S+). Using eqs. ((35l [36)l . it is easy to see that this is 
the line ABS+u = (1 + Sc^ - 9S+)(S - S+). Using A = 8/27, S = 1 and = 1/3 gives = 0.975913 in this linear 
approximation. The numerical continuation of the straight line is easy, since any initial condition close to P_|_ will 
quickly settle down to U. Such a construction using the parameter values above gives Sk = 1.02545, showing that 
the slope of W decreases marginally beyond the linear analysis. 

Solving for the distant part of U is equivalent to finding the physical fiows at long times. For this, we examine 
u ^ 1. From the phase space analysis we see that S decreases as u increases when both are in the physical region. 
Hence, in the second of eqs. we can neglect the term in with respect to the others. Also, the term in (1 -l-Cs)S 
can be neglected with respect to the remaining two terms. Then the system of coupled equations to be solved is — 



^^[l)-[^'^A-B:i"^)- (3^) 

Initial conditions uq and So arc imposed at = in order to match this asymptotic solution with the full solution. 
Unfortunately, even this equation is too complicated for an analytic solution. 

We do not need all solutions of cq. (j37p . as it happens. Since </? < 4/3 for hydrodynamics to apply, we must have 
S < 4/3x = A. For the Boltzmann fluid, then, S < 8/27. Of course, this does not restrict all physical trajectories to 
approach U from below; trajectories can approach this curve from above, but they correspond to a different B. An 
explicit asymptotic form is easy to write down when S <C 3. Then S can be neglected in the equation for u, and one 
gets— 

(3~c^)/4 

To 



[S + 3-cl] /4 



u{t) = ""of") implying e(T) = eq 
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S(r) = cxp 



4:Bu{t) 



Since these forms are asymptotically valid, the constants Sq, uq and eq are free parameters which ensure that the 
asymptotic solution matches the exact solution at large r. The forms above are not to be extrapolated down to small 
r. In this limit one recovers Bjorken scaling, i.e., the boost-invariant solution of the ideal gas. Not only is this a 
satisfactory mathematical result, it could also be of physical relevance, if late freezeout occurs. For later reference we 
note that at late times one has 

using the asymptotic expansion of the exponential integral. The solution illustrated in Figure [8] shows that, as a 
result of viscosity, the expansion is slower than Bjorken, so that the energy density is diluted less rapidly. As a result, 
the initial energy density, as inferred from an observed final energy density, is very much smaller than the Bjorken 
estimate. 

From the phase plane analysis, it is clear that physical flows with initial conditions lying below Sjj are attracted 
to U from below. For such solutions S increases initially before decreasing. Such solutions have been exhibited in 
[Tll.[T^. Other initial conditions for physical flows give rise to monotonically decreasing solutions for S. The solution 
exhibited in Figure [5] is the critical solution, U, which separates these two types of solutions. The long time behaviour 
of all solutions is arbitrarily close to U. 



3. The ELNS approximation 

As discussed before, the ELNS limit of the equations can be obtained by dropping the term in the derivative of ny, 
and then introducing an arbitrary scale of time called Ttt- In this approximation, the equations for the Boltzmann 
fluid become 

U0 = ^u [S + 3~ cl] , and BuS = A. (40) 

Substituting the expression for S obtained from the second equation into the first, the equation can be easily integrated 
with the initial condition u{6 = 0) = 1, to give 

u{e)=e'('-^/'^m, 'M^e-^'f\9), where /(^?) = 1 + ^^^ri^ (l - e^^/^-i)^) , (41) 

and B — 1 + c^. This is in the form of the Bjorken solution modified by a factor which goes to a constant at large 6. 
As before, the ELNS solution matches the IS solution at times much larger than the intrinsic time scale r^. 



4. Entropy production 



For a Boltzmann fluid the entropy density is proportional to a 
into the form 



u^exp(— 36'). One can then manipulate eq. 



da 



3 



(42) 



where we have used the value = 1/3. Since S decreases with time, it is clear that at sufficiently late time the factor 
S — 4/3 becomes negative, and hence the entropy density must decrease. The total entropy scales as E = crexp^, 
since the spatial volume element picks up a scale factor of r from the metric. For this quantity we find the equation 



de 4 ' 



(43) 



which is positive definite since the factors on the right hand side are all positive. Hence the total entropy must 
increase. Using the expression for S in eq. ([55]) . one finds that InS can be written in terms of Meijer-G functions. 
Using instead the asymptotic expansion of S in eq. (j39p . one obtains the simpler long-time limit 

M \ /ToX2/3' 

In contrast to the simple fluid, where the entropy increases without bound, the entropy of an expanding Boltzmann 
fluid reaches a flnite limit. The solution is shown in Figure [SI 



I;(t) = S(to) exp 



(44) 
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FIG. 9: The normalized entropy E as a function of r/r,r. Note the saturation at large r. The numerical solution is for A = 8/27 
and S = 1. The upper curve is for S"(0) = 8/27 and the lower for S°(0) = 2/27. 



B. Sound waves 




FIG. 10: (Color online) Solutions of the linearized equations for fluctuations around the scaling solution in a Boltzmann fluid. 
The normalized temperature, T, the quantity, S, and y are shown for k — 1/4 (red), 1/2 (green), 1 (blue), 2 (purple) and 3 
(orange) for S°(0) = 8/27. The other initial conditions are u^{0) = 1/V2, S^(0) = y^{0) = 1/2. 



The linearized equations for fluctuations in a Boltzmann fluid are 

i(3-c? + SV-^(S-S>X + J-°S\ 

-BS%^ - [(S0)2 - 58° - A]yl - {Bu'' + 28° - B)S\ (45) 

where we have made the expansion f{9, rf) — f'^{9) + Af^{0, rf) where / is any of u, y and S. As before, — 0, and 
vP and S" are obtained as the solution of eqs. ()34|) . and we can Fourier transform in rj to examine the evolution of 
each mode, k. Explicit solutions for and S" can be written only in the long-time limit, when these tend to the 
unstable manifold lA. However, some physically interesting phenomena are likely to occur before this time. Hence, 
the method of choice is to numerically solve eqs. p4l [45]) together for the five variables at the leading and first order 
in A. The initial conditions for and S° have been discussed before. Those for u^, y^ and can be chosen to lie 
between —1 and 1 in order for A to give the right order of the magnitude of fluctuations. 

The results obtained using — 1/3, A — 8/27, B = 1 are shown in Figure [TUl For large values of k, there are 
quasi-periodic solutions. For smaller values of k the solutions are overdamped. The critical value, feg, which separates 
damped and oscillatory solutions depends on 8*^(0). The solutions do not change qualitatively if B is changed by an 
order of magnitude. From the solutions displayed, it is clear that a damping exponent can be extracted. 
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FIG. 11: (Color online) The power spectrum of fluctuations after evolution from an initially flat spectrum at 6 — 0.8 (red), 1 
(green), 1.25 (blue) and 1.5 (black). The numerical solution is obtained for = 1/3, B = 1, A = 8/27 for four different values 
of S''(0). The other initial conditions are (0) = -0.95, y^(0) = -0.5 and S^(0) = 0.3. Note the difference in the scales. 

A quantity which encapsulates the physics, and is perhaps better suited to making a connection with experiments 
is the power spectrum. In Figure [TT] we display the evolution of the power spectrum starting from an initially flat 
spectrum. Qualitatively, the behaviour is reminiscent of the simple fluid examined earlier. Evolution produces peaks 
in the power spectrum. The positions of these peaks evolve with time — moving to smaller k due to the redshifting 
discussed earlier. The position and magnitude of the peaks depend very strongly on initial conditions and Cg- For the 
relation between the power spectrum studied here and the correlation function in [2^, see Section III. 

The connection with ELNS hydrodynamics is made, as before, by dropping the term in the time derivative of 
TTv- As we have mentioned before, this results in the hyperbolic equations turning into a parabolic set. As a result, 
fluctuations are diffusively damped, and do not turn into propagating sound waves. This is seen in numerical solutions, 
and will be dealt with more completely in the next section. 

VI. A CONFORMAL FLUID 

As discussed in Section II. A, a conformal fluid differs from a Boltzmann fluid only in the T-independent value 
X ~ eTjr/riv = Svra, where a — Tt^T. The hydrodynamics of such a conformal fluid can then be taken over from that 
of the Boltzmann fluid with the simple replacement A — 4/3% — A/Qna. In a strongly interacting fluid one might 
expect a — t^^T to be small, and hence A to be large. The separation into boost-invariant and fluctuation equations 
and the analysis of each, is very similar to the details presented in Section V. For this reason, we do not repeat the 
material here, but only point out the differences. 

The boost-invariant flows follow eqs. As discussed in Section V.A, the global structure of the flows is determined 
by the three flxed points P^,, P_ and P+. When A is large enough, the roles of the fixed points P_ and are 
interchanged. However, physical flows are governed by the unstable manifold of P+ and its behaviour remains 
unchanged. In particular, one may take over the expressions for the long-time behaviour of the unstable manifold. 

The boost-invariant flows, starting from the same initial conditions, as x changes, are shown in Figure 1121 The 
late-time behaviour of the trajectories in ip'^ = ^'iid is independent of the value of A. This is clear by using 
eqs. (159]) to write 
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FIG. 12: (Color online) Boost-invariant flow for conformal fluids in the plane of ip^ and starting from the same initial 
conditions, with B = 1, but with varying X) i-S-, -4. As A increases from 1 to 9 in steps of 1 {A increases from top to bottom), 
the long-time behaviour remains universal, but is approached faster. 

This universality is a consequence of the structural stability of the hydrodynamic equations for this class of fluids. 
From the flgure one can also see that for ^ « 1, (p^ increases before decreasing. Thermalization can be said to occur 
only when the solution enters the physical domain < 4/3 for the last time. 




20 40 60 80 100 20 40 60 80 100 

^\ ^\ 

FIG. 13: (Color online) Exact numerical solutions for (a) the energy density and (b) total entropy along the unstable manifold 
for conformal fluids. In (a) the line in black shows the exact solution for ^ = 8, the line in gold is a Bjorken solution matched 
to it asymptotically, and the line in blue is the corresponding solution in ELNS hydrodynamics matched asymptotically. In (b) 
the coloured lines show the time-evolution of E as ^ varies from 1 to 9 in steps of 1, the larger the A, the higher the curve. 
The line in black is the prediction from the ELNS solution matched to the energy density. 

From eq. it is clear that for all conformal fluids the amount of entropy generated during the flow, proportional 
to S, has a finite upper bound. However, this bound increases exponentially with A, starting from the initial value 
5^(''o) ~ independent of A. Note that this means that at fixed value of rjv/s, the maximum entropy production has 
exponential dependence on the inverse relaxation time, l/rjr. 

Since the longitudinal fiow can be continued to times of order tt , and tt is independent of any material property 
other than c^, for conformal fluids, the ratio tt/t,^ could become large as Ttt decreases. One might then expect that 
ELNS evolution should suffice to describe the system. Figure [T^ illustrates several interesting points. First, the late 
time evolution of the energy density can be approximated by an appropriately tuned Bjorken solution, the tuning 
parameter being the initial energy density. However, as shown in the figure, this gives a gross over-estimate of the 
energy density; in the example, by more than two orders of magnitude. Second, a solution of ELNS hydrodynamics 
can be tuned to reproduce the energy density at late times. Once this tuning is performed, it also reproduces the 
total entropy production at late times, and hence furnishes a far superior description of the fiow than the Bjorken 
solution could. Note, however, that the ELNS solution has larger entropy production than the true IS solution at 
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FIG. 14: (Color online) The power spectrum for the energy density, at 9 = 0.8 (red), 1.05 (green), 1.25 (blue) and 1.5 (black), 
starting from a constant unit distribution with initial conditions S''(0) — min( 1,A), u^(0) = -0.95, y^iO) = -0.55 and 
S^(0) = 0.30 for aU k. The successive figures are for (a) A^2, (h) A= 1, (c) A ^1/2, (d) A ^1/3. 



initial times and that there are quantitative lacunae in this approximation even at t = AOttt- With lower A, we have 
seen that these discrepancies are larger. Finally, we note that the maximum entropy production occurs in the very 
early stages of the flow, and probes of this stage of the expansion would best discriminate between different values of 
r^. A different approach to extracting is advocated in [l7|. 

A possible discriminant between ELNS and IS hydrodynamics is the fate of fluctuations around the scaling solution. 
An analysis of fluctuations can be performed numerically using eqs. (|45p . as before. At small k all solutions are 
overdamped, as can be seen by investigating the fc — s- limit, as before. At large k the fluctuations develop into 
damped propagating waves. The evolution of the power spectrum of the fluctuations in energy density, starting from 
a uniform spectrum is shown in Figure [Til For generic initial conditions, increasing A, i.e., decreasing r^r, seems to 
damp fluctuations faster. However, the equations become stiff for A > 2 and the numerical solutions are hard to 
extract for the interesting case of iV = 4 SYM theory which yields A = 8.69. 

However, in that case, we can take another approach. The asymptotic solutions exhibited in eqs. ()381 I39|) are 
reliable for large 6. One can investigate the fate of fluctuations around the scaling solution at late times by inserting 
the asymptotic formulae into eqs. (|45p. The asymptotic solutions can be written as = mq exp(p6') and S° = 
{A/ Buq) exp(—pd), where p = 1 — B/A. We expand eqs. (|45)l in powers of 2; = exp(p6'), and retain all non-negative 
powers of z in the equations. This gives 



where 




M = ifcMi +Mo, 



Mo 




4 ^ 



-uoBeP^ 



(47) 



Note that the set and can be solved independently of u^, and this last variable is then driven by the others. 

For orientation, let us examine some analytical approximations first. As fc — > 0, one may set M = Mq. The 
equations are then exactly solvable, and yield 



S\t) = Sjexp 
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FIG. 15: (Color online) The late-time power spectrum for the energy density, at 9 = 0.8 (red), 1.05 (green), 1.25 (blue) and 1.5 
(black), starting from a constant unit distribution in (a) IS hydrodynamics with ^ = 8 and initial conditions ii^(O) = —0.95, 
y^{0) — 0.95 and S^(0) = —0.75 for all k and (b) ELNS hydrodynamics with the same initial values of and . 

At large k, neglecting Mq, one finds that the solutions for and are oscillatory with frequencies ito = ±iky/ A/ B. 
The solution for is entirely driven by yi, and hence is oscillatory with a growing amplitude. The growth exponent 
precisely matches the growth exponent of hence /u^ is purely oscillatory in this approximation. Damping arises 
with the inclusion oi Mq. Then the equations are not autonomous, and one could hesitate to treat the matrices Mi and 
Mq as time independent. However, by choosing k to be large enough, one may be able to treat exp(p0) as a constant 
over many periods of oscillation. The simplest approach to computing the damping exponent then is to treat the 
problem in perturbation theory. This is best done, as before, by assuming that ioj = ikXi + Aq + 0(\/k), introducing 
this expansion into the characteristic equation for M and then solving for Ao,i. As expected, Ai = -il^/{A/B) and 0. 
The damping exponent is — Aq. We find that y^ and have the common damping exponent [1 — + uoe^^]/2. The 
temperature fiuctuations, /vP are also damped. 

The equations of ELNS hydrodynamics can be recovered from those of IS hydrodynamics by the method which 
we have explained at length earlier. In this case, the ELNS equations are obtained by first setting S° — A/Bu^ 
(see eq. [57]) . an approximation which is equivalent to the late-time solution in eq. and reducing the last of 

eqs. (j47|) through the further identification = Ayl^/B. Then the equation for y^ turns into a diffusion equation, 
which has a completely different character from three equations of IS hydrodynamics which together give rise to wave 
propagation. The qualitative difference between the two kinds of hydrodynamics is very clear from the evolution of 
the power spectrum of temperature fluctuations, as shown in Figure [TSl The propagation of damped sound waves in 
IS hydrodynamics is signalled by the formation of beats, and its absence in the case of ELNS signals the diffusive 
nature of the dynamics. 

VII. CONCLUSIONS 

In this paper we examined freely expanding ffuids using causal viscous hydrodynamics [19] in the longitudinal 
1-dimensional approximation. We chose coordinates appropriate to a description of longitudinal flow, performed a 
tensor analysis and reduced the tensor equations of hydrodynamics to three scalar equations. The hydrodynamic 
modes are described by three scalars, the energy density, e, shear viscous part of the energy-momentum tensor, Try, 
and a parametrization of the velocity, y. Details of this procedure are given in Appendices A and B. 

Material properties expected of the QCD fluid were discussed in Section 11. A, and three models of the fluid were 
put forward. We gave the name "simple fluid" to a model in which Cs, rjv and r^r are temperature independent. Other 
toy models, which go by the names of Boltzmann fluid and conformal fluids, have = 1/3 and x = ^Tt^/Vv, both 
independent of the temperature. Different such fluids are distinguished by the value of x- For the same value of 
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rjvl s, one can have any value for x\ large values of the latter corresponding to gaseous fluids and small values to more 
liquid- like behaviour. 

In Section II. B we performed an analysis of the symmetries of the hydrodynamic equations and extracted dimension- 
less scaling variables which are combinations of the hydrodynamic variables and material properties. These express 
laws of physical similarity, and can be related, in the non-relativistic limit, to the familiar dimensionless variables, 
i.e., the Reynolds number, Re, the Mach number, M, and the Knudsen number, K . For Boltzmann and conformal 
fluids, the condition that x does not depend on the temperature implies the combination KYle/M is constant. We 
discussed how IS hydrodynamics provides a self-consistent description of thermalization. 

In scaling flow of an ideal fluid, the initial energy density entirely goes into driving the expansion. This leads to the 
Bjorken expansion with its characteristic power law behaviour, e oc l/r^ (where B = c^). In contrast, in a viscous 
fluid, some part of this energy density is dissipated, leading to entropy production. Viscous effects are sufficiently 
strong in the simple fluid that the power law is modified to e ~ l/r asymptotically (section IV. A). As a result, the 
total entropy of the system increases without bound (section IV. A. 2); we find S (x r^/'* asymptotically. 

We made a phase plane analysis of flows for conformal fluids, including the Boltzmann fluid (see Section V.A.I). 
Our analysis shows that the long-time behaviour of any physical flow is arbitrarily close to a special solution of the 
hydrodynamic equations which we call the unstable manifold, lA. We investigated this solution in detail (see Sections 
V.A.2 and V.A.4), in particular, the approach to the Bjorken solution and rates of entropy production. 

In the Boltzmann fluid, Bjorken scaling is recovered asymptotically (Section V.A.2). However, the initially slower 
expansion means that the energy density is diluted slower; Bjorken estimates of the initial energy density, given the 
final, are therefore over-estimates. The total entropy is saturated reasonably early (section V.A.4). Such behaviour 
also holds for conformal fluids (section VI). However, for fixed rjy/s, the saturation value for the entropy depends 
exponentially on the relaxation time r^r. The rate of growth of the entropy does not depend uniquely on rjv / s, since 
the parameter t-^ also plays a role. 

When the relaxation time for the shear stress, Ttt is very small in units of 1/T (i.e., the fluid behaves like a liquid), 
then the late time solution for a conformal fluid is reasonably well approximated in ELNS hydrodynamics by tuning 
a free parameter. This free parameter is equivalent to the initial energy density (section VI). After performing such 
a fit, the entropy density can also be fitted provided an appropriate unit of time is chosen. Hence, the extraction of 
initial conditions and material properties from observed final data remain as equivalently hard problems in ELNS and 
IS hydrodynamics. Furthermore, the initial conditions inferred from a given final energy density in the two kinds of 
dynamics differ by a large factor. Hence, in making inferences about the system produced in heavy-ion collisions, one 
must be careful to use the right type of hydrodynamics. 

For an ideal fluid, the fluctuations around the scaling solution are overdamped for spatial Fourier modes k < kg = 
(1 — c1)/2cs- For larger k, the fluctuations become damped sound waves. The energy density in the fluctuations is 
damped as a power law in r (section HI). These qualitative features persist in IS hydrodynamics of viscous fluids 
(sections IV. B, V.B and VI). When k is small enough, the solutions are overdamped. At larger k damped sound waves 
are obtained. The scaling solutions are therefore stable against small fluctuations. The power spectra of fluctuations 
contain interesting structure, which, if observable, could give information on fluid properties or initial conditions. 

Although the scaling solutions obtained using IS hydrodynamics can be well approximated at asymptotically late 
times by the solutions of ELNS hydrodynamics (sections IV. A. 1, V.A.3, VI), fluctuations behave completely differently 
in these two kinds of dynamics. In the latter, there are no sound modes: all fluctuations are diffusively damped. In 
IS dynamics, however, modes with large enough k give rise to true sound waves. These two kinds of behaviour are 
easy enough to distinguish through power spectra of the energy density (Figures [TH [15]). Similar effects are also to be 
seen in the fluctuations of y. Note the rather precise analogy of acoustically produced peaks in the power spectrum 
of the temperature fluctuations in heavy-ion collisions with those in the blackbody radiation in the early universe. 

One object of great modern interest in heavy-ion collisions has been the coupling between hard and soft particles, 
for example jet-quenching. An effective theory description of this is to treat the soft particles through fluid dynamics 
and describe the hard particle as forcing a shock wave in this medium. The behaviour of acoustics is an important 
ingredient in such an analysis. In view of the difference between ELNS and IS hydrodynamics in their treatment of 
fluctuations, they could potentially give rise to different predictions for such observables. 

The purpose of this paper was to explore the kind of phenomena that arise in IS causal viscous hydrodynamics, 
and find qualitative features which are different from either ELNS hydrodynamics or ideal gases. The changes in the 
scaling solution are features which distinguish IS hydrodynamics from ideal gases. The propagation of fluctuations 
around these solutions as sound waves gives a qualitative distinction between IS and ELNS hydrodynamics, since 
these modes are diffusively damped in the latter case, and do not propagate. 

Comparison with data, and fits to initial conditions and material properties are outside the scope of this paper. 
These questions require a proper treatment of radial and anisotropic flows, as well as the incorporation of bulk viscosity 
(due to its importance near and below Tc). These questions are left to the future. 
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APPENDIX A: TENSOR DECOMPOSITIONS 
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TABLE I: Contractions of the symmetric basis tensors. The entries are the contractions of the tensors in the row and column. 
The contraction in the first column is the trace. 



In longitudinal flow there are only two vectors intrinsic to the problem — the timelike u and the spacelike v. One 
could construct two more spacelike vectors to complete a tetrad, but since these vectors are completely arbitrary, no 
physics can depend on them. In order to express the rank-two tensors which enter the hydrodynamic equations, we 
can only use combinations of u and v, and the metric tensor. 

The symmetric rank- two tensors can be chosen to be the projections L^" = u'^u'^, V'^'^ — — u^w", and Yj^^ — 
gtif „ ^A"' _ ytJ'i' In addition there is the traceless symmetric tensor A^'' — u^v'^ -f u'^v^^ , and the antisymmetric tensor 
A^'^ = u'^v'^ — u'^v^. The double contraction of the antisymmetric and any of the four symmetric tensors vanishes. 
For notational convenience we also define the spacelike projection A^" — g^" — L^'' = S'^'^ -I- V'^'^ . Contractions of this 
with the other tensors can be worked out using Table IH Any rank-2 symmetric tensor, which arises in consideration 
of longitudinal flow, can be expressed as a linear combination of L, V and E (or alternatively, L, V and A). Any 
similar rank-2 antisymmetric tensor can only be a scalar multiple of A. 

The vorticity tensor is defined as 

ojt^- ^A'^^A'^^^idpu^-d^up), (Al) 

and is spacelike (orthogonal to u) and antisymmetric by construction. Due to its antisymmetry, u^'^ has to be 
proportional to A. However, A is not orthogonal to u, so the only possible constant of proportionality is zero. In 
other words, vorticity vanishes for longitudinal flow. Another way to understand this is to note that w^"^, being an 
antisymmetric spacelike tensor, is equivalent (technically, dual) to an axial vector which is spacelike. Such a vector 
can be constructed by the three-dimensional vector product of two vectors. For longitudinal flow, there is only one 
spacelike vector v which is intrinsic to the problem. Hence one cannot construct an axial vector. A tedious proof of 
the vanishing of the vorticity can also be given by direct manipulation of the definition. 

In the presence of viscous shear but vanishing bulk viscosity, one of the hydrodynamic variables is the dissipative 
part of stress tensor, tt^''. The fact that it is symmetric can be derived from the symmetry of the stress tensor. Since 
it expresses shear, it is orthogonal to u. It is traceless since we have assumed bulk viscosity to vanish. Hence one can 
write uniquely. 



2 



(A2) 



i.e., for longitudinal flow, the shear tensor is equivalent to a single scalar. Note that the tensor multiplying Try has 
the following properties — 



Iy.'"' ^ Iv^"" - ^A^"" and [V" - l^^niVnu ^ l^n^] = 
2 2 2 ^ 2 '"^ 2 2 



We will need to consider derivatives of the basis tensors. This is facilitated by considering first the derivatives of the 
unit vectors. The condition u'^m^ = 1 yields u^^Du^ — 0, and therefore the orthogonality of u and v. Now D^u^ can 
be written as a linear combination of u and v. It is easy to work out that D^u^ = S^u^ + {DS)v^ . Furthermore, one 
can show that Dv^^ — Su'^. Two other special cases of interest for longitudinal flow are d^v^^ = S and v'^Du^ = —9. 
For later use we also set down the actions of £> = u^'d^t and D = w^dp on any scalar field / — 

Df = fr cosh 2/ + -f,, sinhy and Df = fr sinhy + -f,, cosh?/. (A4) 

T T 

Using the orthogonality of u and v, we find the derivative 



V^u'" = d^M^ = e, where = A^d^. (A5) 
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Another derivative that appears repeatedly in the hydrodynamic equations is- 



(A6) 



where the notation (■ ■ •) denotes the traceless, symmetric part. Since this tensor is also orthogonal to u, one should 
be able to write (see eq. IA2|) 



V ~h ' 



where T) is to be determined. Multiplying through by y^i,, one immediately finds 



(V^u,) = -2 ( v^Dux + I 



T/ V 



Now we are ready to write the derivatives of the basis tensors. 
Since covariant derivatives of the metric tensor vanish, one has 

dxA^^u -{d\u^t)ui, - {d\Uu)u^,. 

We examine its projections parallel and orthogonal to u. One of the parallel projections is 

DA^u = —SA^u. 

One of the orthogonal projections is 

A^VpA'"" = Ai^^A^dxA^P = 0. 

For V^jy one has 

dxVi_,^ = -{dxVi_,)vi, - {dxVu)v^. 
A parallel projection which we will use later is 

DVf_,i, = -{Dv^)v^ - {Dv^)vf_, = ~SA^„. 
One of the orthogonal projections that we need is 

Ai^WpVP = [S- dpvP]v''. 



(A7) 
(A8) 

(A9) 
(AlO) 
(All) 
(A12) 
(A13) 
(A14) 



Since Vpv'^ — —1, one has v'^Dvp — 0, i.e., Dv^ is parallel to u. As a result, one finds Dv^ = {u^Dvx)u^, from which 
the last form of the derivative follows. 



APPENDIX B: THE EQUATIONS OF LONGITUDINAL HYDRODYNAMICS 

The dynamical equations are supplemented by the equation of state, which provides a relation between e and p, and 
hence determines the entropy density s = (e + p)/T. Since the the hydrodynamic equations are valid only for a fluid 
in local thermodynamic equilibrium, or so close to it that linear response theory works, one may use the equation of 
state to eliminate one of e and p from hydrodynamics. A toy equation of state that we shall use is 



p(T)=c2e(r), 



(Bl) 



i.e., the speed of sound is independent of the temperature. The only value of which is strictly temperature 
independent is = 1/3. This is the appropriate value to use when the bulk viscosity has been neglected, since both 
are consequences of conformal symmetry. Since we use temperature independent in this paper, we use the above 
value whenever numerical work is performed. 
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1. The equation for energy 

The equation for energy in [TJ] is 

De = -(e+p)V^^,A< + i7r^'^(V^u,). (B2) 

The identity in eq. (|A3|) . used along with eq. (|A8p gives 

TT^'iV^u,) = 27rve. (B3) 
Then, using the relation in eq. (IA5|) . and the equation of state, one can write 

De = - [(1 + cl)e - Try] 9. (B4) 

2. The momentum-balance equation 

The general form of the momentum-balance equation given in is 

(e + p)!??/^ = V> - A^VpTT'"'' + Ti>"'Du„. (B5) 

Note that each term is orthogonal to u. Using the definition of v and the decomposition of the shear tensor in eq. 
\2\\. we find that 



tt'^'^Du^ = Stt^'^'v^ = Snvv^- (B6) 
Using the derivatives of the basis tensors in eqs. (jAlll IA14p . we can write 

A^Vp^^" - A^A^(dAvry)[U'^'' - ^E^"] - ^-^vv^'ld^v^ - S]. (B7) 

The last term drops out because d\v^ — S. Putting all this together, we can reduce the tensor equation to 

S{€ + p)v'' ^V^'p ~ {dxT:v)[V^^ + Sttvv" . (B8) 

Contracting with any spacelike tensor apart from v would yield only terms in the directional derivative along that 
vector. But for longitudinal flow each such derivative is separately zero. Thus, the only nontrivial equation is obtained 
by contracting the above equation with v. This gives the entropy equation for longitudinal flow, 

S'(e + p-7ry) + 5(p-7ry) = 0. (B9) 

3. The equation for the shear tensor 

In the equation for the shear tensor is given as 

t^A^A^Dtt"'' + ^^-^ = ??^(V^w") - 2r,^"(^u;j;), (BIO) 

where r^r is a relaxation time, rjv is the coefficient of shear viscosity, — A'"'c?^, and w"'^ is the vorticity tensor. 
Since the vorticity vanishes for longitudinal flow, we drop the last term. 
Using eq. (jASp the term in the viscous coefficient becomes 

w(v'^0 = ^??.e[y'^"^-isn. (Bii) 

A^A^Dtt"^ = (Li7ry)[y^'' - ^EH, (B12) 

since the projectors are orthogonal to the derivatives of the basis tensors in eqs. (jAlOl [5T3| . The equation for shear 
then reduces to the scalar equation 

r^(i?^y) + Try = ^T^v-e. (B13) 



For the first term, one can write 
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APPENDIX C: EVOLUTION OF THE POWER SPECTRUM 



For the equations of sound we assume that a Fourier transformation decouples the individual Fourier modes, 
each mode, the ODEs to be solved can be written in the form 



For 



(CI) 

where x is a vector and M a matrix. For ideal fluids, eqs. (|15|) can be written as a system of two coupled autonomous 
equations, so that Af is a 2 x 2 matrix. In the other cases (eqs. [301 SSI and its reduction in the ELNS limit) the 
equations are not autonomous, i.e., the time variable appears explicitly in M. However, the solution of eq. (|C1|) is 
straightforward and involves a matrix exponential. The asymptotic behaviour of the solution is controlled by the 
eigenvalue with the largest real part. We will assume that all the eigenvalues of M have negative real part (at all 
times, in the non-autonomous case), so that at long times the solution decays. We are interested in whether transients 
grow. 

Consider the real positive quantity P{0) = x^Ax where ^ is a fixed matrix independent of 9. Let A4 = M'^ A + AM . 
Then 



where 



(C2) 



in terms of the normalized eigenvectors, v^, and eigenvalues. A;, of A^. Since P is real, A is Hermitean. As a result, Ai 
is also Hermitean, and its eigenvalues are real. Transient growth in P can take place for some initial conditions if and 
only if M has at least one positive eigenvalue. The largest growth in P occurs when x is parallel to the eigenvector 
of A4 with largest (positive) eigenvalue. 







FIG. 16: (Color online) The power spectrum of fluctuations in energy density for the ideal fluid as a function of k and 6, for 
/3 = and a = 0, n/4, Tr/2 and 37r/4 (from left to right). Note the lack of transient growth for fc = 0. 






FIG. 17: (Color online) The power spectrum of fluctuations in energy density for the simple fluid in ELNS hydrodynamics as 
a function of k and 9, for f3 — and a = 0, n/4, n/2 and 37r/4 (from left to right). Note the lack of transient growth for k — 0. 

When A is chosen to be a multiple of the identity, then the eigenvalues of A4 are multiples of the real parts of the 
eigenvalues of M . Then, since the multiplying constant is positive, and we have assumed that the eigenvalues of M 
have negative real parts, there can be no growth in P. In all other cases, transient growth of P is possible. 

In our application, we choose A to be a projection operator on the first component of x. As a result, AM is 
just the top row of M with other rows set to zero, and A4 is obtained by Hermitizing this. Now the eigenvalues 
of Ad have no simple relation with those of M. One can show that in general the rank of A4 is two, i.e., there are 
two non-vanishing eigenvalues. Unless all the off-diagonal terms in M are zero, one of the generically non-vanishing 
eigenvalues is positive and the other is negative. As a result, there is always transient growth in P. This is the reason 
for the peaks in Figures [HI [H and [TT] 



26 



We complete the analysis of transient growth of P for the ideal fluid (see eas. ITS)). In this case we have 

~ *0 ^ ^ ' with eigenvalues A = ±kB. (C3) 

Transient growth can take place. We can parametrize all initial conditions by an angle a and a phase /3, by choosing 
X = (sinQ;exp[i/3],cosa). The evolution of the power spectrum of the energy obtained with specific initial conditions 
from numerical solutions of eqs. (|15p are shown in Figure 1161 

For the evolution of the power spectrum of in a simple fluid, one has 



3 „ \ ... , . „ 1 




tVI = I 4ifc j , with eigenvalues A = 0, ^ 



-4± V25 + 144fc2 



(C4) 



where we have taken — 1/3. Since one of the eigenvalues is positive, transient growth occurs. Note that the mode 
fc = can also display transient growth. This analysis can be adapted to that for the Boltzmann fluid by replacing 
the constants in M by appropriate functions of 9. However, the conclusions regarding transient growth carry over to 
that case. 

For the power spectrum of in the ELNS approximation to the simple fluid, one flnds that 
2B \ 1 



M = -\ i_ 2 \ ^ with eigenvalues \ = B± —^/Qc^^B^ + 16^(1 - c^,f. (C5) 

lik^ J ■iCs 

Interestingly, the diffusive term, c^k^, drops out of the transient analysis for the energy density (it does appears in the 
transient analysis for y). For fc = transients do not grow. This seems to be the main distinction between transient 
growth of energy density in ELNS and IS dynamics. 

The numerical results in Figure [17] seem similar to the ideal fluid case at first sight. However, the diffusive character 
of the underlying equations is manifested in two ways. First, at fixed fc, if one observes the time evolution, then one 
sees only a single instance of transient growth, unlike the quasiperiodic behaviour of in the ideal fluid. Second, at 
sufficiently large fc there is no transient growth, unlike the ideal fluid. 

In the late time limit for a Boltzmann or conformal fluid one has, from cq. (|47p . 



2p 



4 
( 

^eP^ 




M=\'-^^eP'> ' with eigenvalues A = 0,p± - Jl6p2 + (i + B2p)y2g2pe_ (eg) 



4 

Since one of the eigenvalues is positive, transient growth can take place. 



[1] A. Hosoya, M. Sakagami and M. Takao, Ann. Phys., 154, 299, (1984); 

A. Hosoya and K. Kajantie, Nud. Phys., B, 250, 666, (1985); 

G. Baym et al, Nucl. Phys., A, 525, 415C, (1991); 

S. Jeon and L. G. Yaffe, Phys. Rev., D, 53, 5799, (1996); 

P. Arnold, G. D. Moore and L. G. Yaffe, J. H. E. P., 11, 001, (2000); 

G. Aarts and J. M. M. Resco, J. H. E. P., 4, 053, (2002); 

P. Kovtun, D. T. Son and A. O. Starinets, J. H. E. P., 310, 064, (2003). 
[2] C. Eckart, Phys. Rev., 58, 919, (1940). 

[3] L. D. Landau and E. M. Lifshitz, Fluid Mechanics, Elsevier, New Delhi (2005). 
[4] J. D. Bjorken, Phys. Rev., D, 27, 140, (1983). 
[5] G. Baym et al, Nucl. Phys., A, 407, 541, (1983). 

[6] S. Kagiyama, A. Nakamura and A. Minaka, Prog. Theor. Phys., 76, 171 (1986). 
[7] D. H. Rischke, S. Bernard and J. A. Maruhn, Nucl. Phys., A, 595, 346, (1995); 

P. F. Kolb, J. SoUfrank and U. W. Heinz, Phys. Rev. C , 62, 054909, (2000); 

D. Teaney, J. Lauret and E. V. Shuryak, Phys. Rev. Lett. , 86, 4783, (2001); 

P. Huovinen et al, Phys. Lett., B, 503, 58, (2001); 

T. Hirano and K. Tsuda, Phys. Rev. C , 66, 054905, (2002); 

P. F. Kolb and R. Rapp, Phys. Rev. C , 67, 044903, (2003). 
[8] See for example, the treatment in P. J. E. Peebles and Bharat Ratra, Rev. Mod. Phys. , 75, 559, (2003). 



27 



[9] P. Danielewicz and M. Gyulassy, Phys. Rev., D, 31, 53, (1985); 
D. Teaney, Phys. Rev. C , 68, 034913, (2003). 
[10] M. Prakash, M. Prakash, R. Venugopalan and G. Welke, Phys. Rep., 227, 321, (1993). 
[11] A. Muronga, Phys. Rev. Lett. , 88, 062302, (2002) and Phys. Rev. C , 69, 034903, (2004). 
[12] A. Muronga and D. H. Rischke, nucl-th/0407114 

[13] U. Heinz, H. Song and A. K. Chaudhuri, Phys. Rev. C , 73, 034904, (2006). 

[14] R. Baier, P. Romatschke and U. A. Wiedemann, Phys. Rev. C , 73, 064903 (2006), and JVud. Phys., A, 782, 313 (2007). 

[15] M. Lubinsky and E. Shuryak, arXiv:0704.164Zl 

[16] P. Romatschke and U. Romatschke, arXiv:0706.1522| 
[17] A. Dumitru, E. Molnar and Y. Nara, arXiv:0706. 22031 
[18] H. D. Weymann, Am. J. Phys., 35, 488, (1967); 

D. C. Kelly, Am. J. Phys., 36, 585, (1968); 

N. G. van Kampen, Physica, 46, 315, (1970). 
[19] W. Israel, Ann. Phys., 100, 310, (1976); 

J. M. Stewart, Proc. Roy Soc, A 357, 59, (1977); 

W. Israel and J. M. Stewart, Ann. Phys., 118, 341, (1979); 

W. A. Hiscock and L. Lindblom, Ann. Phys., 151, 4 66, (1983). 
[20] T. Koide, G. S. Denicol, Ph. Mota and T. Kodama, |hep-ph/0609l"l7 
[21] M. Ichiyanagi, Prog. Theor. Phys., 84, 810, (1990); 

T. Koide, Phys. Rev., E, 75, 060103(R) (2007). 
[22] A. Muronga, Phys. Rev. C , 76, 014909, (2007). 
[23] A. Muronga, Phys. Rev. C , 76, 014910, (2007). 

[24] L. D. Landau and E. M. Lifshitz, Classical Theory of Fields, Elsevier, New Delhi (2005). 
[25] H. Kouno, M. Maruyama, F. Takagi and K. Saito, Phys. Rev., D, 41, 2903, (1990). 
[26] R. Baier and P. Romatschke, nucl-th/0610108 

[27] R. Jackiw in Current Algebra and Anomalies, eds. S. B. Treiman et al, p. 168, World Scientific, Singapore, 1985. 

[28] S. Gupta, Pramana, 61, 877, (2003). 

[29] B. Beinlich et al, Eur. Phys. J., C 6, 133, (1999); 

Y. Namekawa et al, Phys. Rev., D, 64, 074507, (2001); 

A. Ali Khan et al, Phys. Rev., D, 64, 074510, (2001). 
[30] P. Arnold, C. Dogan and G. D. Moore, Phys. Rev., D, 74, 085021, (2006). 
[31] R. V. Gavai, S. Gupta and S. Mukherjee, Phys. Rev., D, 71, 074013, (2005). 
[32] D. Kharzeev and K. Tuchin, arXiv:0705.4280, 
[33] M. P. Heller and R. A. Janik, hep-th/07032^ 
[34] S. Gupta, Phys. Lett, B, 597, 57, (20 04); 

H. B. Meyer, eprint arXi v:0704.180T1 [hep-lat]. 
[35] S. Gupta, hep-ph/0507210: 

R. S. Bhalerao et al, Phys. Lett., B, 627, 49 (2005). 
[36] I. S. Gradshteyn and I. M. Ryzhik, Tables of Integrals, Series, and Products, Academic Press, San Diego, 2000. 
[37] http : //mathworld. wolf ram. com/Mei jerG-function.html 

[38] V. I. Arnold, Ordinary Differential Equations, The MIT Press, Cambridge, Massachusetts, 1973. 

[39] See, for example, A. C. King, J. Billingham and S. R. Otto, Differential Equations, Linear, Nonlinear, Ordinary, Partial, 
Cambridge University Press, Cambridge UK, 2003. 



